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We study Landau damping in dilute Bose-Einstein condensed gases in both spherical and prolate 
ellipsoidal harmonic traps. We solve the Bogoliubov equations for the mode spectrum in both of 
these cases, and calculate the damping by summing over transitions between excited quasiparticle 
states. The results for the spherical case are compared to those obtained in the Hartree-Fock 
approximation, where the excitations take on a single-particle character, and excellent agreement 
between the two approaches is found. We have also taken the semiclassical limit of the Hartree- 
Fock approximation and obtain a novel expression for the Landau damping rate involving the time 
. dependent self-diffusion function of the thermal cloud. As a final approach, we study the decay of a 

condensate mode by making use of dynamical simulations in which both the condensate and thermal 
cloud are evolved explicitly as a function of time. A detailed comparison of all these methods over 
a wide range of sample sizes and trap geometries is presented. 
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I. INTRODUCTION 



CO One of the more challenging problems in the study of Bose-Einstein condensation (BEC) in trapped atomic gases 

concerns the finite temperature dynamics of a condensate interacting with a noncondensed, thermal component. Of 
particular interest in this regard is the study of collective oscillations of the condensate, where one would expect to 
see the influence of the thermal cloud on both the frequency and damping of the modes. Experiments which address 
| these properties have been performed in a number of laboratories around the world 0, 0111,^1, an d provide an ideal 
£^ . test-bed for the development of theories of Bose-condensed gases at finite temperatures. 

One theoretical approach that can be used in the analysis of these experiments is the ZNG theory [jj. In this 
CD , approach, the condensate is described by means of a generalized Gross-Pitaevskii equation while the dynamics of 
the thermal cloud is determined by a Boltzmann-like kinetic equation. The full numerical implementation of this 
7-H ! theory has recently been presented and its use in the analysis of a variety of experiments has appeared in a 
J> ■ series of papers H, 0, El • An important part of this previous work involved a detailed study of the collisional 
dynamics which is responsible for the equilibration of the system with respect to both energy and particle number. 
In this paper, however, we shall focus on an aspect for which collisions are of secondary importance, namely Landau 
damping. This damping mechanism is dominant in the low density systems studied experimentally 0, U 3)- In 
this so-called collisionless regime, mean-field interactions mediate the transfer of energy from the condensate to the 
' thermal component, leading to the damping of condensate collective modes. 

The subject of Landau damping in dilute BECs has been explored by several authors. Early work in an uniform 
p ■ gas was conducted by Hohenberg and Martin llll for low temperatures, while results for higher temperatures were 
obtained by Szepfalusy and Kondor 12|. Liu |l3j used results obtained for the uniform gas to estimate the damping 
in the case of a trapped gas while Pitaevskii and Stringari 0] and Fedichev et al. developed expressions for the 
Landau damping rate that explicitly accounted for the trap geometry. Similar results were derived by Giorgini [Tri Il7| 
who also considered the Baliaev damping process whereby a collective mode decays into two lower energy excitations. 
The latter can be ignored, however, for the low-lying collective modes in a trapped gas because of energy conservation 
£j , restrictions. Semiclassical results for Landau damping were also obtained using a dielectric formalism by Reidl et al. 
fl8| , who found quite good agreement with the results of experiment [| . The ZNG theory also provides a description 
of Landau damping and applications of the theory to a number of different condensate collective modes 7, 8, 9] found 
damping rates that agreed well with experiment. 

The perturbation theory approach of Refs. 0, ITrl Il7| provides an expression for the Landau damping rate that 
involves a sum over transitions between pairs of Bogoliubov excitations. The first fully quantum mechanical evaluation 
of this expression was performed by Guilleumas and Pitaevskii (which we shall hence forward denote as G-P I) [l9| . 
for a spherically-symmetric trapping potential. We should also mention that a related calculation was performed by 
Das and Bergeman pof. obtaining similar results to G-P I. In the ZNG theory, the description of Landau damping is 
quite different since it involves the dynamical evolution of the condensate in the presence of a thermal cloud treated 
semiclassically at the level of the Hartree-Fock approximation. The relation of this approach to that used in G-P I is 
far from evident, although it is clear that the physical basis of the two is the same. To investigate this question, we 
performed numerical simulations in an earlier paper [|| for the system studied in G-P I and in fact found quantitative 
agreement. One of the aims of the present paper is to compare the two methods over a much larger range of condensate 
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sizes, which means that we have had to repeat the kind of calculations performed in G-P I. By doing so, we are then 
able to investigate the consequences of using Hartree-Fock, as opposed to Bogoliubov, excitations in the Landau 
damping calculation. A reformulation in terms of Hartree-Fock excitations has the added benefit of allowing the 
semiclassical limit to be taken in a straightforward way, thereby completing the connection to the ZNG theory. 

Perhaps an even more interesting question concerns the damping rate in anisotropic traps. This issue arose in 
experiments on the transverse breathing mode in a highly elongated trapped gas Q, where the damping rate was 
found to be anomalously low - around an order of magnitude smaller than one would expect in similar circumstances 
on the basis of both experiment and theory. In Ref. |9( we simulated this experiment and found that the low damping 
rate could be explained as a consequence of an in-phase collective oscillation of both the condensate and thermal 
cloud. On the other hand, if the thermal cloud is not set into oscillation, the Landau damping rate is found to have a 
value typical of that seen in other experiments. This result apparently contradicts the conclusion of Guillcumas and 
Pitaevksii (G-P II) [2l| who performed a Landau damping calculation for an infinitely long cylindrical condensate. 
They found a damping rate which was extremely low, around an order of magnitude smaller than the experimental 
results, and two orders of magnitude smaller than our simulation results for a stationary thermal cloud. A natural 
question concerns the validity of modelling the experimental geometry by an infinite cylindrical trap which changes 
the nature of the excitation spectrum and presumably the damping rate. In this paper we address this question by 
repeating the G-P I calculations as a function of the trap anisotropy and compare these results with those obtained 
from our simulations. Our tentative conclusion is that the damping rate has an anomalous dependence on anisotropy 
in the limit of the infinite cylinder, so while the G-P II result is correct, it does not apply at the experimental 
anisotropy. 



II. THEORY 



A. Landau damping of Bogoliubov excitations 



For purposes of completeness it is useful to summarize the Bogoliubov theory and the way in which a condensate 
collective mode is damped by thermal excitations. In a second-quantized notation, the Hamiltonian is given by 

H = J drft(r) + Wr)) &r) + § j dv^(r)^(r)^(v)^(v) (1) 

where 14 x t(r) is the trapping potential and g — A-Kah 2 /m is the interaction parameter defined in terms of the s-wave 
scattering length a. If the field operator is expressed as 

V> = $o + ^, (2) 

where $o is the condensate wavefunction and ifi is the excitation (or fluctuation) field operator, the Bogoliubov 
approximation is generated by expanding H to second order in ijj. The terms linear in if) vanish if <!>o satisfies the 
Gross-Pitaevskii equation 

-h 2 V 2 \ 

+ V CKt + gn c0 $0 = , (3) 



2m 

where n C Q = |<I>o| 2 is the condensate density. At this level of approximation, the condensate does not interact directly 
with the thermal excitations as determined in the Bogoliubov theory. 

The resulting Bogoliubov Hamiltonian is then diagonalized by means of the Bogoliubov transformation 

4>(r) = J2(u i (r)a i -v*(r)ai) . (4) 

i 

The quasiparticle amplitudes Ui and Vi satisfy the Bogoliubov equations 

Lu t - gn c0 Vi = EiUi 

Lvi - gricoUt = -EiVi , (5) 

where we introduce the operator L = — h 2, V 2 /2m+V ex t+2gn c o — /z. The orthonormality of the quasiparticle amplitudes 
is specified by the relation 



dr [u* (r) Uj (r) - v* (r) Vj (r)] = 5 l] . (6) 
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Apart from a constant, the Bogoliubov Hamiltonian is given by 

h b = Y. E A a i > ( 7 ) 

i 

where the excitation energies Ei correspond to condensate collective modes. In thermal equilibrium, these modes are 
populated according to a Bose distribution fo(Ei) — [exp((3Ei) — l]^ 1 with (3 = (fe^T) -1 . 
The cubic terms in the expansion of H, 

V {3) = g J dr $ (ftrprp + ft flip} , (8) 

couple the Bogoliubov excitations and lead to a mechanism for their decay. In particular, a mode which has been 
excited, and therefore no longer in equilibrium with the other thermal excitations, will experience a decay known as 
Landau damping. To represent this situation, one defines a state in which the mode occupation n osc is large compared 
to the equilibrium value. The energy in this mode, E osc — huj osc n osc , then decays as n osc goes to equilibrium. The 
transition rate from the initial nonequilibrium state to any other state can be determined by perturbation theory [hH 
and the average rate of change of the energy in this mode is found to be 

iP = "IT E \ A v\ 2 (h - fM E i ~ E * ~ ^so) , (9) 

y 

where the transition matrix element is 

Aij = 2g J dr $ [( u iU* ~ viu* + v.v*) u osc - (uiU* - UiV* + v t v*) v osc ] . (10) 

For future reference we note that the Hartree-Fock approximation is recovered by setting Vi (but not v osc ) in 110|) 
equal to zero and determining Uj from the first of the Bogoliubov equations in JSJ with Uj = 0. 

The damping rate T is defined according to 2r = — E osc /E osc (which implies that the amplitude of the mode decays 
as e" rt ). Thus, 



r = \ 2 {h - fi)6{Ej -Ei- hco osc ) . (11) 

ij 



This expression is the basis of the calculation of Landau damping as carried out by Guilleumas and Pitaevskii for both 
spherically symmetric and cylindrical traps |l9l l2lj | . In order to display the results of our calculations it is convenient 
to follow their notation 

r 



where 



y^^ijS(uJij -CJqbc), (12) 



~ta = Ki—\M 2 {h-fil (13) 



is a "damping strength" associated with each possible transition, with frequency difference = (Ej — Ei)/K. 

The calculation of the damping then consists of solving the Bogoliubov equations (jSJ for it osc , t> osc and lu osc for the 
classical mode of interest, together with the corresponding quantities for the thermally-populated excitations. The 
matrix elements can then be evaluated for transitions between each pair of excitations, which in turn yields 7^ . 
More details of these calculations will be given in the next section, while the calculation of the total damping rate T, 
including how we deal with the delta function in (|12f) . will be presented in Section III. 



B. Anisotropic traps 

In performing the above calculations, a significant saving in computational effort can be realised by considering 
a spherically symmetric trap. Not only does (|1(J|I effectively reduce to a one-dimensional integral, but the spherical 
symmetry also reduces the number of excitations that have to be considered. In particular, if we label the modes with 
the usual quantum numbers {n,l,m}, then states corresponding to a given / are (21 + l)-fold degenerate, and Aij 
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need only be calculated explicitly for one of these states. This spherical case was first studied in Ref. [19j, and we will 
repeat some of these calculations in Section 3A in order to compare various approximations. However, here we are also 
interested in the damping in anisotropic traps which is the most common experimental situation. It turns out that 
calculations for anisotropic traps, even those with axial symmetry, are much more involved since the degeneracy of 
excitations with different m quantum numbers is lifted. As a result, determining the Bogoliubov excitation spectrum 
requires the diagonalization of much larger matrices than in the spherical case. The corresponding increase in the 
number of possible transitions, along with the fact that evaluation of the Ay matrix elements now involves integrating 
over both radial and angular variables, significantly increases computational overheads. In the remaining part of this 
section we shall describe how we tackle these problems by making use of a spherical harmonic basis. The calculation 
for an isotropic trap is then a simplified version of this general scheme. 

As a first step, we find the ground state condensate wavefunction as detailed in Ref. j2^. This makes use of 
a set of normalised basis functions ipnimi?) = Rni(r)Yi m (8, </>), which are eigenfunctions of the Hamiltonian ho = 
— ft 2 V 2 /2m + Vo(r). The potential Vo(r) is the I = component of the total axially-symmetric effective potential 
V(r) = K x t(r) + gn c (r) when expanded in terms of Legendre polynomials, V(r, 6) = Vi (r)Pi (cos 8). We take the 
external harmonic potential to be of the form T4 x t(r) = (muj 2 ^/2)(x 2 + y 2 + X 2 z 2 ), where A is the anisotropy parameter. 
By defining a nonspherical perturbation Al^(r) = V(r)— Vo(r), and expanding an arbitrary solution of the GP equation 
as 4>(r) — ^2 n ia n iip n im{^), one can set up a matrix equation for the coefficients a n i- The condensate wavefunction 
is then given by the lowest-energy even-parity solution &o(r) — VAT c o (r), which in turn can be used to calculate 
the condensate density n c (r). Updating the effective potential V(r) then leads to a new matrix equation, which is 
solved to yield an updated wavefunction. The calculation is repeated until $ converges to the correct condensate 
wavefunction, with the lowest eigenvalue giving the chemical potential fi. In the process, one also generates a full set 
of eigenfunctions Q (r) for the Hamiltonian h = — h 2 V 2 /2m + V(r) — fi with eigenvalues e a . 

Once we have this information, the Bogoliubov equations can be solved by introducing the expansion "0 t + ( r ) — 
X) Q c "Va( r ): where ipf(v) = Ui(r) ± Vi(r). The index i labels the different excitations and includes the azimuthal 
quantum number m. Again, this leads to a matrix equation, which can be solved to obtain the eigenfunctions ^^(r) 
and energy eigenvalue Ei for each mode. For the purposes of the Landau damping calculation, which we move onto 
next, it is convenient to express these functions as 

^+(r) = J2 d niRm(r)Y lm (e^), (14) 
V>"(r) = Y. e niRm{r)Y lm (0^), (15) 

nl 

by making use of the expansions of (j) a in terms of the ip n im- 

For the particular mode of interest (i.e., the one for which we are calculating the damping) we can define the 
mode densities, 5n + (r) = &o(r)ip+ sc (r), and (5n - (r) = <I>o(r)-0~ c (r). Using this notation, the Landau damping matrix 
elements become 

At] = f Jdr [Sn-(^+*^ + Wj*1>7) + 5n + ^+^~ - </>">+)] . (16) 

In evaluating this integral, it is convenient to expand the mode densities (where the mode has a particular m-number, 
which we denote m) in terms of spherical harmonics 



^±(r) = £^^*nf(r)y fe (M), (17) 



where 



21 + 1 



2k 



^0 



Sn n r ) = \l-^-l W I dBBmBY^faffln^T). (18) 



Substituting lfH)l. (JTHJl and |T7|) into (JT^J gives 

drr 2 £ Rnl^Rnn^l^d^+Se^e^B^r) 



Al > ~ 2 



nln'V 



+ (d$e%-e$d%)B+ ij (r)], (19) 
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where 



/ 2V -4-1 

4 y » - ^^^^{Tl'mmijl'mm^m^nfirl (20) 

in which ilil^mi'm^ll^'ms) is the Clebsch-Gordan coefficient [2^ |. and where mj = m, + m. The calculation of the 
Landau damping then involves evaluation of (|19|l for each pair of excitations with a difference of energies in a particular 
interval i? m i n < Ej — Ei < E max . In this paper we shall focus on the Landau damping of modes with m = and even 
parity, so that selection rules connect excitations with Am = and the same parity. Unlike the spherically-symmetric 
case, however, excitations with different m values must be summed explicitly. Nevertheless, since the contributions 
from the ±mj excitations are identical, a two- fold saving can be gained for rrii =/= 0. Our results for the breathing 
mode as a function of anisotropy will be presented in Section 3B. 



C. Hartree-Fock Approximation 

In this section we obtain an expression for Landau damping in the Hartree-Fock (HF) approximation. Although 
this can be done trivially by setting — in (|ll)f) . it is useful to rederive the result using an alternative method for 
a number of reasons. First it clearly shows that Landau damping is associated with the work done on the thermal 
cloud by the dynamic mean field of the oscillating condensate. Second, this reformulation allows one to take the 
semiclassical limit in a straightforward way, thereby facilitating a comparison with the scmiclassical simulations that 
we have performed on the basis of the ZNG theory. 

Within the HF approximation, the condensate and thermal cloud are treated as two distinct components. The 
condensate evolves dynamically according to a time-dependent GP equation while the thermal cloud responds to the 
condensate mean field. This picture is quite distinct from the Bogoliubov approach in which the 'condensate mode' is 
distinguished from the 'thermal excitations' only through the assumed different occupations of the respective modes. 
Common to both pictures, however, is the absence of a thermal cloud mean field acting back on the condensate. Such 
effects arc included in the ZNG theory and are crucial for the satisfaction of the generalized Kohn theorem |5j • 

For small deviations from equilibrium, the time-dependent condensate wavefunction is given by 

*(r,t)=# (r) + *$(r,t), (21) 

where the fluctuation S$ is obtained from the linearized GP equation. The fluctuating part of the condensate density 
is then 

6n c (r, t) = $ (r) [<S$(r, t) + <5$*(r, t)} (22) 

which gives rise to a dynamic mean field 2gSn c (r, t) acting on the thermal cloud. Thus the thermal cloud experiences 
a perturbation 

H'(t) = J dr2gSn c {r,t)n(r) , (23) 

where n(r) represents the density operator of the thermal component. The linearized GP equation has solutions of 
the form 

5*(r, t) = Uosc (r)e— _ < sc (r)e^ 4 , (24) 

where the u and v are, apart from a normalization, the same Bogoliubov amplitudes introduced earlier and w osc is 
the frequency of the particular condensate mode of interest. In terms of this wavefunction, the condensate density 
fluctuation is given by 

5n c (v,t) = $ (r)fc c (r) e -'^'+c.c. 

= £n-(r)e- iuwt + c.c. (25) 

where 



V'oscM = ?w(r) - tw(r) 



(26) 
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With this form of the condensate density fluctuation, the perturbation given in (|23|) has a harmonic time dependence 
and can be treated by means of time-dependent perturbation theory. One thus finds that the time-averaged rate of 
change of the thermal cloud energy is given to lowest order in the perturbation by the expression 

2u osc (2g) 2 [ dr [ dr'6n-*(r)x"{r,r',Lj osc )5n-(r') (27) 



dt 

where x"(r, r', ui) is the imaginary part of the time Fourier transform of the thermal cloud density response function 

X (r, r', t-t')= l -6{t - t')([n(r, t),n(r', i')])o . (28) 

The angular brackets, (...)o, denote an equilibrium expectation value. Treating the thermal atoms as independent HF 
excitations, the spectral density is given by the expression 

X "(r,r» =nJ24>*(r)h(r)$(r')ct> i (r')[fi ~ fjWi " * ~ M ( 2 9) 

ij 

where (f>i(r) is the HF wavefunction (as obtained from JSJ by setting Vi = 0) for the i-th state with energy e^, and fi 
is the thermal Bose occupation number. Substituting this expression into (|27|l . we find 

— = 2nLU osc ^ \ A ij \ 2 (fi - fj)&{ £ 3 - Si- ftWosc) (30) 

ij 

with 

A l3 = 2. 9 / dv & (r)cf>* (r)5n- (r) . (31) 



As stated earlier, this matrix element can be obtained from ({101) by simply setting Vi = 0, but retaining v osc which is 
needed to properly define the condensate density fluctuation. 

The rate of change of the condensate energy E osc is of course the negative of (|30|l . To extract a damping rate 
we must divide this by the energy of the mode which is determined by the amplitude of the condensate density 
fluctuation. If u and v are normalised according to JfJJl, this energy is just huj osc . 



D. Semiclassical HF 



Our purpose in this section is to show how the quantal HF damping can be reduced to a semiclassical form. This 
will allow us to make contact with the semiclassical ZNG simulations. The response function in l|28() can be written 
X(r, r', t-t') = 2i0{t - t')(j)(r, r',t- t') with 

0(r, r', t -t') = ^r([n(r, t), n(r', f)]) . (32) 
This correlation function can be expressed in the form [24| 

1 r p - 

<p(r,r',t) = - d\(h(r,t-ih\)h(r')) a (33) 

2 Jo 

which allows the semiclassically limit for the dynamics to be taken straightforwardly. Setting h — > 0, we have 

0(r,r',t) = ^(n(r,t)n(r'))o (34) 
where now n(r, t) is a classical density variable. The Fourier transform of (|34|l gives 

cj>(r,r',cj) = ^ / ^e Mt (n(r,t)n(r')) . (35) 
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To evaluate the correlation function in we consider a thermal cloud containing N atoms. The density of the 
cloud at time t is thus given by 

N 

f»(r,t)=2*(r-r < (t)) (36) 

i=l 

where rj(i) is the position of the i-th particle at time t. Thus, 

(n(r, t)n(r')) = *£{6(r - r t (t))5(r> - Tj (0))) . (37) 

Since the atoms in the cloud are independent and equivalent, we have 

(n(r,t)n(r')) =G s (r,r',t)+ho(v)h (r') (38) 
where we have defined the self-diffusion function [2{| 

G s (T,v',t) =iV<5(r-r 1 (t))5(r / -r 1 (0)))o (39) 

and no(r) is the equilibrium thermal cloud density. Since the product of equilibrium densities is time-independent, 
it will not contribute to the Fourier transform in l|35|) at finite to and can be dropped. Combining these results, the 
semiclassical (sc) Landau damping rate is given by 

— = W^r^o) (40) 

where f(oj) is the Fourier transform of the function 

f(t) = { dr ( dr'Sn-*(r)G s {r,r',t)Sn-(r') . (41) 



We thus see that the quantity determining the damping rate is an appropriately weighted spatial average of the 
self-diffusion function G s . The latter can be evaluated analytically for some simple model situations which provide 
insight into its behaviour. However in the present situation, the thermal atoms are moving in the combined potential 
of the trap and the mean field of the condensate and G s can only be determined numerically. It is then easier to deal 
with f(t) directly. Using the definition of G s , f(t) can be expressed as 

f(t)=N{5n-*(r 1 (t))8n-(r 1 (0))) , (42) 

which is convenient for numerical evaluation. 

We note that ri(t) is the position of an atom at time t starting from some initial position ri(0) = ro with momentum 
Pi(0) = po- The initial phase point is distributed according to the equilibrium Bose distribution /o( r Oj Po) 5 normalized 
to N, and so the correlation function takes the form 

f(t) = |^P°/o(ro,po)<5n-*(r 1 (ro ! po;t))<5n-(r ) (43) 

Thus, f(t) can be determined by averaging over all possible trajectories starting from some initial phase point. In 
practice this can be done by considering an ensemble of Nt test particles distributed in phase space according to 
/o(i"o,po) and then performing the discrete sum 

N Nt 

Nt is chosen sufficiently large to ensure that statistical fluctuations in /(i) are acceptably small. Examples of such 
calculations will be given later. Finally, a numerical Fourier transform of f(t) at to = lu osc yields the semi-classical 
damping rate. 

The value of /(0) can also be evaluated directly as 

/(0) = J drn (r)|fc-(r)| 2 . (45) 

This serves as a check of the sum over test particles described above. In addition, we note that f(t) tends to a finite 
limiting value for t — > oo. This limiting value can be subtracted from f(t) in the evaluation of /(w) since this only 
affects f(ui) at u> = 0. 
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E. ZNG Simulations 

Our final method for evaluating the Landau damping is based on A^-body simulations |6( in which collisions are 
neglected. In this numerical scheme, the condensate wavefunction is evolved in time by numerically solving the 
generalized GP equation using a split-operator FFT method, while the evolution of the thermal cloud is calculated 
by sampling with a large number of test particles, which are evolved classically. This is in fact equivalent to solving a 
collisionless Boltzmann kinetic equation for the thermal cloud, where the thermal excitations are described according 
to the semiclassical Hartree-Fock (HF) approximation (i.e., as single particle excitations in a quasi-uniform gas 0). 
The dynamic mean field of the thermal cloud then gives rise to decay of the condensate oscillation, which is identified 
with Landau damping. Full details of the numerical methods are given in Ref. [(|. 

The first step in the calculation is the self-consistent determination of the condensate wavefunction together with 
the thermal cloud density. The simulation is then initiated by exciting the condensate in such a way as to realise the 
relevant mode of study. It is straightforward to show that this can be achieved by adding ip~ BC (r), as obtained from 
the solution of the Bogoliubov equations JSJ), to the condensate ground state wavefunction. The subsequent evolution 
of the mode of interest can then be monitored by evaluating the moments (x 2 ), (y 2 ), and (z 2 ) of the condensate 
as a function of time. These are calculated numerically using (x)(t) = J dr x\&{r,t)\ 2 ■ A plot against time yields 
a damped oscillation, which is fit to a sinusoidal function with an exponentially decaying factor exp(— Tt). Each 
simulation extends for a time of uj±t =10, where the relatively short time is chosen so as to avoid driving the thermal 
cloud too far out of equilibrium. This point is discussed more fully in 6J. The damping rate, T, is then extracted from 
the simulation data, which can then be compared to the results of our other calculations. It should be emphasized 
that within this approach the damping of the condensate oscillation is a direct result of the dynamic mean field that 
the thermal cloud exerts on the condensate. In other words, no damping would occur if the condensate were to simply 
oscillate in the static equilibrium mean field of the thermal cloud. 

III. RESULTS 
A. Isotropic Traps 

We now move on to describing the results of our calculations, where we begin with the case of the Bogoliubov 
method for an isotropic trap (A = 1, u>o = oj± = 2tt x 187 Hz). The specific system being considered is a gas of Rb 
atoms with a scattering length a — 5.82 x 10 -9 m. The calculation follows closely that of Guilleumas and Pitaevskii 
(G-P) [l^. The results are conveniently displayed as a histogram which plots the damping strengths, 7^ , versus 
the excitation frequency = (Ej — Ei)/h. Fig. [Ja) shows such a plot for N c — 2.5 x 10 5 condensate atoms at a 
temperature of ksT j \x — 1.5, where the chemical potential takes its Thomas-Fermi value /1 = 29.9?kJo- The spectrum 
features several high peaks, which correspond to transitions between low- lying excitations with large overlaps (i.e., 
large Aij matrix elements) together with large population factors, fi — fj. Following G-P, we shall refer to these peaks 
as "resonances" . They are not expected to play a significant role in the Landau damping process, unless one happens 
to lie very close to the mode frequency uj osc . For this reason we shall ignore these transitions when calculating the 
Landau damping, and instead focus on the small amplitude quasi-continuous "background". 

A conceptual difficulty arises in the calculation in that the excitation spectrum in l|ll|) consists of discrete delta 
functions. Taken literally, this would imply that the Landau damping rate would either be zero or infinite depending 
on the location of the oscillation frequency ui osc . Following this difficulty is overcome by replacing each delta 
function in the background spectrum by a Lorentzian A/{(27r)[(cJij — oj OS c) 2 + A 2 /4]}. Physically, this accounts 
for the fact that the excitations will themselves have a finite lifetime in a more rigorous theory. As a result, the 
Landau damping rate is determined by averaging over many peaks weighted according to their proximity to the mode 
frequency. The width factor, A, however is somewhat arbitrary, and as shown in Fig. [21 the result for T will vary 
with A. Nevertheless, one can see that the variation is weak when A/u>q lies between 0.05 and 0.20. We fit the data 
in this range to a straight line, and extrapolate back to A = to yield an estimate of the Landau damping rate. The 
result as a function of temperature is shown in Fig. [31 for different numbers of condensate atoms. One sees in the 
plots the typical rapid increase of the rate at low temperatures followed by a nearly linear temperature dependence 
at the higher temperatures. Our results for N c = 5 x 10 4 are in complete agreement with those of Guilleumas and 
Pitaevskii [l9j. 

We have repeated the calculation for HF thermal excitations following the method discussed in Sec. Ill CI The 
spectral density in this case is shown in Fig. ^b), for the same parameters as in (a). The main difference between 
the two plots is the absence of the strong resonances seen in the Bogoliubov calculation. However, the background 
spectrum is remarkably similar, especially in the vicinity of the mode frequency lu Q sc- Insight into why this is happening 
is provided by comparing the Bogoliubov and HF excitation frequencies, as plotted in Fig. 01 for N c = 5 x 10 4 . As 
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1.8 2.0 2.2 2.4 2.6 

FIG. 1: Histogram showing results of Landau damping calculations for N c = 2.5 x 10 5 87 Rb atoms at a temperature of 
kBT/n = 1.5. The height of the bars represents the damping strength -yij of each transition, with frequency uiij plotted on the 
horizontal axis, (a) shows data for the Bogoliubov modes, while (b) is the corresponding HF calculation. Note that in (a) the 
highest peaks extend far beyond the vertical range, since we are most interested in the "background" spectrum. The mode 
oscillation frequency for this radial breathing mode is uj osc = 2.235tJo. 
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FIG. 2: Variation of the damping rate (in units of the mode frequency u osc ) with the Lorentzian width A (in units of the 
trap frequency), for A = 1 and N c = 5 x 10 4 . Each curve represents a different temperature, with (moving from bottom to 
top) ksT/fi = 0.4, 0.8, 1.2, 1.6, 2.0 and 2.4. The straight lines are fits to the data between A/u>o = 0.05 and 0.20, where the 
intercept with the j/-axis gives the estimated damping rate. 



found previously |26( , the two spectra are completely different at low energies and angular momenta, but converge at 
high E and I, demonstrating that the excitations take on a single-particle character j2j|2^ in this region. It is these 
excitations that are responsible for the majority of the background transitions in Fig. ^ which explains the similarity 
between the HF (calculated in the same way as described above) and Bogoliubov damping rates plotted in Fig-El The 
excellent agreement with the Bogoliubov results for a wide range of condensate sizes is in contrast with the uniform 
Bose gas where one finds a significant difference between the two approaches. The reason for this is that the 
Bogoliubov spectrum for a uniform gas only approaches the single-particle HF form at temperatures much larger than 
the chemical potential, while in a trapped gas surface excitations with high multipolaritics arc also important, even 
at low temperatures. We thus arrive at the important conclusion that the 'collective nature' of the excitations has 
little bearing on the determination of Landau damping in trapped Bose gases. 

It is also of interest to compare the Bogoliubov results to those of our semiclassical simulations based on the ZNG 
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FIG. 3: Damping rate (in units of u osc ) plotted versus temperature, fc_gT//i, where \i is the Thomas-Fermi chemical potential 
for each number of condensate atoms, N c . Results from a Bogoliubov calculation are plotted as solid lines, while those for the 
HF approximation are plotted with dashed lines. Results for N c = 10 4 , 5 x 10 4 , and 2.5 x 10 5 are shown. 
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I 

FIG. 4: Excitation spectrum for Bogoliubov (lines) and HF excitations (bullets), for a spherical condensate with N c = 5 x 10 4 . 
For each mode the energy E and multipolarity I are plotted. The horizontal line at pi = 15.8/kJo represents the chemical 
potential of the condensate. 

theory (as described in Section 2E). We plot the damping rates for both calculations as a function of the number 
of condensate atoms in Fig. [3 for three different temperatures. The behaviour seen for both sets of calculations is 
similar, with a more rapid rate of increase of the damping rate at low N c followed by a less rapid increase at higher 
N c . The data at N c = 5 x 10 4 and N c = 1.5 x 10 5 correspond to those plotted in Fig. 8 of Ref. @, where we 
concluded that the two approaches were in good agreement. The comparison in Fig. [5] shows that the agreement 
persists over a wider range of condensate number, although some systematic differences do exist. The simulation 
rates at low N c tend to be lower than the Bogoliubov damping rate, while at high N c there is a tendency for them 
to be slightly larger. These small differences are presumably due to the different approximations used in the two sets 
of calculations. These include (i) the use of HF as opposed to Bogoliubov excitations in the simulations, (ii) the use 
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FIG. 5: Damping rate (in units of ui OS c) against number of condensate atoms, N c for an isotropic trap. The results from 
a Bogoliubov calculation (open symbols) are compared to those of a semiclassical simulation (solid), for temperatures of 
hsT/n = 1 (circles), ksT/fi — 1.5 (squares), and ksT/fi — 2 (triangles). The dashed lines through the Bogoliubov data serve 
as guides to the eye. Each simulation data point is calculated by taking the mean of damping rates extracted from oscillations 
in the three directions, with the standard deviation giving the error bar. 



of the semiclassical approximation for the thermal cloud dynamics and (iii) the inclusion of the thermal cloud mean 
field in the self-consistent calculation of both the equilibrium and dynamical properties. Concerning (i), we saw in 
Fig. |3| that there is little difference between the Bogoliubov and HF results over a wide range of temperatures when 
the excitations are both treated quantum mechanically. We therefore do not expect the use of Bogoliubov, as opposed 
to HF, quasiparticles to appreciably change the results. 

The effect of the semiclassical approximation itself, however, is less clear. To isolate this effect we make use of the 
formulation outlined in Sec. IID. The calculation is based on Eq. (|40|l with f(t) determined from 144|) by sampling the 
mode density 5n~(r) along classical trajectories. To be specific, we again consider the radial breathing mode in an 
isotropic trap with frequency loq — 2ir x 187 s . The mode density being sampled in this case is shown in Fig. [5] for 
N c = 10 5 . The node in the mode density is of course required in order that the volume integral of the mode density 
vanish. 

We consider first a model in which the thermal atoms move in a purely harmonic potential, that is, we ignore 
the mean field of the condensate. Furthermore, we assume that the phase space density of the atoms is given by a 
Boltzmann distribution. With these assumptions, the self-diffusion function G s (r,r',t) can be evaluated analytically 
with the result p9| 

(3muj 2 (r 2 + r' 2 — 2r • r' cos Lo t) \ 
2 sin 2 ui()t ) 

This function is periodic with period To = 2tt/luq: it recovers its t — value G s (r,r',0) — <5(r — r')rio(r) at the 
times t n = uTq. This is a special property of purely harmonic confinement. When substituted into l|41(l . the resulting 
function f{t) has period To/2 since the mode density Sn~(r) has inversion symmetry. We show in Fig. EJa) f(t) as 
calculated from l|44|) by following the trajectories of an ensemble of test particles at a temperature fc^T '/ ' [i = 1.5. 
The figure confirms the behaviour expected on the basis of the analytic form of the self-diffusion function in i|4rj|) . In 
Fig. CJb) we show a similar calculation of f(t) but now with the condensate mean field included. The distribution of 
thermal atoms is still Boltzmann-like. The strict periodicity seen in Fig. Qa) is washed out as a result of the motion 
of the thermal atoms in the nonharmonic confining potential, Vth(r) = 14 x t(r) + 2gn c (r), although there are clearly 
vestiges of the dominant 2wo frequency of the purely harmonic case. The mean field of the condensate disrupts the 
periodicity of the particle orbits that occurs for harmonic confinement and as a result, f(t) saturates at long times to 
a constant limiting value. The inset of Fig.JJfb) shows f(t) on a smaller time scale; the behaviour seen is consistent 
with f(t) being an even function of time. 

In Fig. 02a) we show f(t) for the same conditions as in Fig. EJb), but with the Boltzmann distribution replaced by 



G s (ry,t)=N(-p*_) exp 
\2n\smujMJ \ 
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FIG. 6: Mode density (arbitrary units) of the radial breathing mode as a function of the radial distance r in units of the 
harmonic oscillator length a^ a — \J h/muio. The condensate contains N c = 10 5 atoms. 
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FIG. 7: (a) The weighted self-diffusion function f(t) for a cloud of thermal atoms in an isotropic harmonic potential with 
frequency loq = 2-k x 187 s" 1 at a temperature of ksT/fi = 1.5. (b) As in (a) but with the mean field potential of N c — 
5 x 10 4 condensate atoms included. The vertical scales are in units of 10 -3 N c /a^ . The Boltzmann distribution used in these 
calculations is normalized to the same number of atoms as obtained with the Bose distribution in Fig. |H| 



the Bose distribution. The behaviour of f(t) is qualitatively very similar; the main difference is the limiting long-time 
value which can be attributed to the different thermal distributions in the two simulations. It is worth commenting 
on what this asymptotic value is due to. If the position of an atom at long times is uncorrelated with its starting 
point, one might expect the self-diffusion function to tend towards no(r)n$(r') / N . This would yield an 'equilibrium' 
value of f(t) of 



f =1 



drSn (r)n (r) 



(47) 
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FIG. 8: (a) As in Fig. □ but with the Bolt zmann distribution replaced by the Bose distribution (solid line). The dashed curve 
shows the effect of including collisions between the thermal atoms. The horizontal asymptote at / eq , Eq. 1471 . is shown, (b) The 
Fourier transform of the weighted self-diffusion function for the two curves in (a): no collisions (solid), with collisions (dashed). 
Both vertical scales have the same units as in Fig. [T] The vertical line indicates the mode frequency at lu osc — 2.235ojo- 



We find, however, that f(t) in our simulations does not tend to this limit. The reason for this is that long-time 
correlations persist since a given atom retains its initial energy in the course of the dynamical evolution. Even if the 
nonharmonic perturbation were to lead to an ergodic distribution on a given energy surface, the equilibrium limiting 
form can only arise if equilibrating collisions between thermal atoms take place. This was confirmed by performing a 
simulation in which collisions are included following the methods discussed in Ref. |(j . The dashed curve in Fig. 02a) 
shows the result of this calculation and we now indeed find that f(t) — > / cq as t — > oo. 

In Fig. IHtb) we show the Fourier transform of f(t). It has peaks at frequencies close to multiples of 2wo a s would 
be expected in view of Fig. OUa). Less obvious is the fact that J(uj) appears to be positive definite. This behaviour, 
however, is necessary since f(uj) is just the semiclassical limit of the spectral density in l|3U|) . It is therefore reassuring 
that this property emerges from our simulations. The dashed curve in Fig.^b) is the corresponding Fourier transform 
when collisions are included. The main difference between the two curves occurs at low frequencies where the collisional 
curve has a Lorentzian peak due to the relaxation of f{t) to / cq on a time scale of to^t ~ 10. 

The semiclassical Landau damping rate is determined by the value of f(u>) at uj — lu osc for the radial breathing 
mode. This frequency is indicated by the vertical line in Fig. E^b) and it is clear that collisions are not important 
in this case. The Landau damping rates calculated on the basis of (|40|l are plotted as a function of temperature in 
Fig.[5]which also displays the Bogoliubov results for comparison. The two sets of calculations are qualitatively similar 
and show a similar dependence on the condensate number N c . However, the semiclassical results underestimate the 
Bogoliubov results at higher temperatures. 

To understand these differences we show in Fig. 1101 excitation spectral densities for the three approximations on a 
much larger frequency range. For the Bogoliubov and HF calculations, the spectral densities are given by (|llfl . with 
lu osc replaced by an arbitrary frequency o>, and convolved with a Lorentzian of width A = O.Iujq. This representation 
is more informative than the histograms in Fig. ^ since it reveals the total spectral weight as a function of frequency. 
The corresponding semiclassical spectral density follows from l|4U|) and is given by 

r S cH = ^c/H. (48) 

It is clear from this figure that the spectral densities in all three approximations are qualitatively similar, with a series 
of broad peaks near multiples of 2ojq. However, both the Bogoliubov and HF calculations have a higher spectral density 
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FIG. 9: A comparison of the Bogoliubov damping rates (lines) with the semiclassical HF approximation (points): N c — 5 x 10 4 , 
solid line and filled points; N c — 1.5 x 10 5 , dashed line and unfilled points. The stars represent the result of a semiclassical 
simulation for N c — 5 x 10 4 as described in the text. 









\ (a) 






1 1 1 






i i i 
1 A (c) 



2 4 6 8 10 

co/co 

FIG. 10: Comparison of the (a) Bogoliubov, (b) HF, and (c) semiclassical HF spectral densities (all having the same arbitrary 
scale) as discussed in the text. The calculations are for a temperature of ksT/fj, — 1.5 and N c = 5 x 10 4 condensate atoms. 
The vertical line at lo = 2.235cjo indicates the position of the radial breathing mode. 



in the range 2ljo < uj < 4wo than the semiclassical result, which accounts for the smaller semiclassical damping rates 
shown in Fig. [!|] The larger quantum spectral densities are associated with low-energy states for which the zero-point 
energy and the effects of barrier penetration are more important. These quantum effects are of course missed in the 
semiclassical HF calculation. However, the overall similarity of the curves in Fig. ^| indicates that the semiclassical 
HF approximation is providing a good qualitative description of the thermal cloud excitations. 

In order to gain more insight into these results, we have performed another simulation which incorporates the same 
physics as the semiclassical HF approximation formulated in Sec. Ill Dl As in the full simulations in Sec. Ill El we evolve 
the condensate using the time-dependent GP equation but do not include the mean field of the thermal cloud acting 
on the condensate. The condensate thus oscillates with a fixed amplitude giving rise to a harmonic perturbation of 
the thermal cloud which itself starts off as an equilibrium distribution and evolves classically in time in the presence 
of the dynamic mean field of the condensate. These are the same ingredients entering the HF perturbation theory 
calculation. We then monitor the thermal cloud energy E t h{t) = Yli [Pi(t)/2 m + ^ifr( r i(*))] > where Vj/i(r) is the 
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equilibrium effective potential acting on the thermal cloud. As a result of the work done by the condensate on the 
thermal cloud, the latter experiences a secular increase in energy which we can identify with the time-averaged rate 
of energy transfer in l|30|) . Some points obtained in this way are shown in Fig. E3 These points agree quite well 
with the Bogoliubov rates, and by the same token, with the full simulations as shown in Fig. |SJ Thus, both types 
of simulations are mutually consistent but deviate from the results based on the semiclassical HF expression in (|40|l . 
Apparently, following the dynamics of the thermal cloud in the presence of an oscillating condensate is not entirely 
equivalent to the evaluation of the equilibrium correlation function in the semiclassical HF expression. At the present 
time we do not have a satisfactory explanation for the observed difference. We speculate that it may be related to 
the nonequilibrium nature of the thermal cloud distribution in the simulations, but we have not been able to find a 
way to check this. 



B. Anisotropic Traps 

We now describe the effect of introducing an anisotropy (A < 1) such that the condensate takes on a prolate geometry. 
The anisotropy lifts the degeneracy with respect to the m quantum number which necessitates the inclusion of a much 
larger number of transitions in our Bogoliubov calculations. As a result, calculations for N c > 10 4 are numerically 
prohibitive. The same is true for large anisotropies, so we are limited here to 0.3 < A < 1. This range, however, is 
sufficiently large to reveal some general trends. 
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FIG. 11: Histogram of Landau damping strength, jij, versus excitation frequency uiij for the m — breathing mode in a 
condensate with N c — 10 4 atoms at a temperature of feT/(i = 1.5. Each plot is for a different trap anisotropy: from top to 
bottom, A = 1, 0.95, and 0.9. The geometric mean of the trap frequency in all cases is LOho = 2n x 187 Hz. 

The effect of imposing a slight anisotropy on the excitation spectrum is shown in Fig. ^2 (for iV c = 10 4 and with 
the geometric mean of the trap frequency Who = w^A 1 / 3 kept fixed at 2n x 187 Hz). One sees that some of the larger 
resonances, corresponding to transitions between the lower-energy excitations, are each split into a series of shorter 
peaks. This is most clearly illustrated for the large peak at uiij/ufho — 1.94, which is associate with a pair of low- lying 
I = 1 modes. The anisotropy splits this peak into three since the m = — 1, and 1 modes are no longer degenerate. 
As the anisotropy increases the peaks continue to separate, as shown by the bottom panel in Fig. ^] for A = 0.9. 

An obvious consequence of this is that the spectrum is now much denser, but one can extract the damping rate in 
much the same way as described above for isotropic traps. The only difference is that we no longer ignore the larger 
peaks (resonances) in our calculation, since these are more numerous and less distinguishable from the background. 
There is therefore no simple criterion that can be used to remove them. The damping rate for each A is plotted in 
Fig. 1121 where one sees a downward trend as the anisotropy increases towards a cigar-shaped geometry (decreasing 
A). The anomalously high damping rate at A = 0.55 is a consequence of a mixing of the radial breathing mode with 
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FIG. 12: Damping rate of the monopole mode (in units of the mode frequency, uiosc) in an axisymmetric trap, as a function of 
the anisotropy, A. The open circles and closed circles represent the Bogoliubov and simulation results respectively. The solid 
line is a straight line fit to the Bogoliubov data, neglecting the outlying data point at A = 0.55. 



another m = mode which happens to possess a similar frequency at this anisotropy. We show the effect of this 
mode crossing on the spectrum by scanning through A = 0.55 in Fig. 1131 Many more high resonances are seen at 
A = 0.55 than to either side which in turn leads to the much higher damping rate shown in Fig. We have here a 
case in which the 'resonant modes' do in fact contribute to the Landau damping. We also believe that the scatter in 
the Bogoliubov data seen in Fig. 1131 can be accounted for in terms of resonant peaks approaching or receding from 
the vicinity of the breathing mode frequency. 
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FIG. 13: Same as Fig. 1111 but for anisotropies of (from top to bottom) A = 0.6, 0.55 and 0.5, showing the origin of the 
anomalously high damping rate at A = 0.55. The shift of the spectrum towards the right is due to the ratio u)±/wh increasing 
as A increases. 

The straight line in Fig.^|is a least-squared fit to the Bogoliubov data (excluding the high-lying point at A = 0.55). 
Extrapolating the straight line fit to A = suggests a damping rate of \/u> osc ~ 0.006 in this limit, which is around 
a factor of 4 smaller than in the spherical case. For comparison, we have also plotted in Fig. ^] the results of 
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semiclassical ZNG simulations. These calculations are again carried out for the same condensate number (N c = 10 4 ) 
and temperature (fe^T//! = 1.5). We note from Fig. [S] that the simulation damping rate at A = 1 lies below 
the Bogoliubov calculation for this low number of condensate atoms. This seems to impart a slightly different A 
dependence as compared to the Bogoliubov results but we still see a downward trend for A < 0.5. In particular, there 
is a rapid decrease in the simulation damping rate for A — > 0. 

These results shed some light on differences found in previously reported calculations 0, I^J. Guilleumas and 
Pitaevskii 21] performed a Bogoliubov calculation for a cylindrical condensate (corresponding to A = 0) and found a 
Landau damping about two order of magnitudes smaller than the typical values in Fig. [5] and an order of magnitude 
smaller than the damping rates observed experimentally for A = 0.065 Q|. On the other hand, we found a 
conventional Landau damping rate from ZNG simulations (as discussed in this paper) that was consistent with Fig. [5] 
but an order of magnitude larger than experiment. The rapid variation of the simulation results in Fig. 1121 for small 
A is a possible resolution of these different theoretical predictions. If the damping rate continues to drop rapidly for 
A — ► 0, one might find a result consistent with that of Guilleumas and Pitaevskii |2l|. Unfortunately, we cannot 
push our calculations to lower A in order to confirm this. It is nevertheless clear from all of our calculations that the 
Landau damping in the cylindrical geometry appears to be an anomalous limit. 

As a final point, we should emphasize that the conventional Landau damping rate discussed here fails to account 
for the recent observations of the damping of the transverse breathing mode in elongated condensates Q . This could 
only be achieved by including the full dynamics of the thermal cloud in the ZNG simulation . Conventional Landau 
damping deals with the oscillation of the condensate in an otherwise equilibrium thermal cloud; the condensate does 
work on the thermal cloud and loses energy to it. However, if the thermal cloud is itself set into motion as is the case 
in the experiments 0|, a change in the damping rate is to be expected. The experimental observations Q can only 
be accounted for if this effect is taken into account 

IV. CONCLUSIONS 

In this paper we have studied Landau damping of Bose-Einstein condensates in isotropic and anisotropic traps and 
have compared a variety of methods for calculating the damping. Calculations based on perturbation theory, involving 
sums over transitions between excited states, show excellent agreement between a Bogoliubov treatment of the exci- 
tations and the Hartree-Fock approximation over a wide range of temperatures and numbers of atoms. These results 
demonstrate that Landau damping in trapped Bose-condensed gases is essentially a 'single-particle' phenomenon in 
that the collective nature of the excitations is unimportant. A semiclassical limit of the HF approximation has also 
been developed, leading to a novel expression for the damping rate in terms of a classical correlation function. A 
detailed comparison with the quantal HF results has shown that the excitation spectra in the two cases are quali- 
tatively similar but differ at a quantitative level. In particular, the semiclassical results tend to underestimate the 
quantum results at higher temperatures. We also explored damping by means of iV-body simulations (as discussed in 
Ref. |3) in two complementary ways. In the first, the condensate wavefunction is evolved using the time-dependent 
GP equation and its interaction with the thermal cloud leads to a decay of the condensate oscillation amplitude. In 
the second, the condensate oscillates with a fixed amplitude and the rate of increase of the thermal cloud energy is 
determined. Both methods give very similar results which agree with the quantum perturbation results over a wide 
range of temperatures and condensate sizes, reinforcing the conclusions of Ref. Q- However, all of these results differ 
somewhat (see Fig.[5J) from the semiclassical HF formulation. At present we have no explanation for these differences. 

We have also studied in detail the damping in cigar-shaped traps using the Bogoliubov approximation and semiclas- 
sical simulations. These results show a decreasing trend in the damping as the trap becomes increasingly anisotropic. 
Over the range of A in both sets of calculations, the damping decreases by about a factor of two, with the semiclassical 
simulations indicating a more rapid decrease for A — > 0. This behaviour may reconcile the differences in the results 
obtained for highly elongated condensates @ and those obtained [2lJ in the cylindrical limit (A = 0). 
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